if (!require('knitr')) install.packages('knitr'); library('knitr')
## Loading required package: knitr
knitr::opts_chunk$set(warning=FALSE, message=FALSE)
# Code by Austin Greene
# Edited by Chris Wall

# Purpose: Generate a map of Waimea Valley, Oahu with elevation data, 
# and site locations which are scaled the mean-annual precipitation 

# Load in packages
require('FedData')
require('latticeExtra')
require('scales')
require('viridis')
require('ggmap')
require('phyloseq')
require('raster')
require('rgdal')
require('RColorBrewer')
require('ggplot2')

# Load in physeq object
physeq1=readRDS("data/physeq1.gz")

# Getting base map via Stamen of location within bounding box. 
# Using Stamen maps which do not require a Google API key
map3 <- get_map(location = c(left = -158.070631, bottom = 21.598356, right = -157.998464, top = 21.655101), zoom=13, source="stamen", maptype = c("terrain-background"))
Waimea_map_3 <- ggmap(map3) # Turn into ggmap graphical object
Waimea_map_3 # Check that it looks good 

# Generating extents of map to pull elevation
bb <- attr(map3, "bb") # Mapp attribute object, from which we extract map extents 
extentB <- polygon_from_extent(raster::extent(bb$ll.lon, bb$ur.lon, bb$ll.lat, bb$ur.lat),
                               proj4string = "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs")

# Confirm bounding box for our elevation data
ggmap(map3) + geom_polygon(data=extentB, fill=NA, aes(x=long,y=lat), color="red", size=3)

# Download elevation data tiles
elev_waimea <- get_ned(template = extentB, label = "ned_waimea", res="1", force.redo = F)

# Make df of elevation data with lat and long included
elev_raster_df <- raster::as.data.frame(elev_waimea, xy=TRUE) #xy True 

# Make a dataframe of geographic location and mean annual precip at each site (from existing metadata)
geo_p=cbind(sample_data(physeq1)$Long, sample_data(physeq1)$Lat, sample_data(physeq1)$rain, sample_data(physeq1)$FieldSite)
geo3 = data.frame(geo_p) #we convert it to a dataframe
colnames(geo3)=c("Long", "Lat","Rain", "Site") # Repair column names
geo3 <- subset(geo3, Site != 9) # Remove Site 9 
geo3$Site <- factor(geo3$Site) # Convert sites to factors, not critical

# Plot elevational data and site locations scaled by precipitation on top of existing ggmap of Waimea Valley
Waimea_map_6_sizescaled <- Waimea_map_3 + # Existing map
  geom_raster(data = elev_raster_df, aes(x=x, y=y, fill=elev_raster_df$ned_waimea_NED_1)) + # Raster elevation data
  geom_point(data = geo3, aes(x=Long, y=Lat, size=Rain), fill="white", color="black", pch=21, alpha=0.5) + # Points for sites, scaled by mean-annual precipitation
  scale_size_area() + # Scale by area. Scaling by size if you delete this line. 
  coord_cartesian() +  # Set to cartesian coordinates
  scale_fill_viridis(option = "plasma", alpha = 0.7) + # Set fill colors to be colorblind-friendly
  coord_fixed(1.3) + # Aspect ratio of cartesian coordinates is 1.3
  labs(x="Longitude", y="Latitude", size="Precipitation (mm/year)", fill="Elevation (m)") + # Labels
  theme_classic() # Classic minimalist theme

Waimea_map_6_sizescaled # Print the plot

# Save the plot as a tiff image with 300 dpi. 
# Setting to 6x6 inches of the plot with the expectation that this is 2x larger than publication size 
# per publication guidelines (expected print size 3x3 in)
ggsave("figures/Map6_SizeScaled.tiff", Waimea_map_6_sizescaled, height = 6, width = 6, units = "in", device = "tiff", dpi=300)

Mantel: fung, bact

# Code by Anthony
# modifed by Chris Wall
# To plot mantel correlation between fungi and Bacteria
require(vegan)
require(ggplot2)

#read normalized Bacteria Phylloseq object
PS_normal=readRDS("data/PS_normal.gz")

#read normalized Fungi Phylloseq object
fungal_PS_normal=readRDS("data/fungal_PS_normal.gz")

#sqrttransform
sqrtfun=transform_sample_counts(fungal_PS_normal, function(x) x^0.5)

#remove the sample missing from Bacteria
sqrtfun=subset_samples(fungal_PS_normal, CollectionID !="WMEA_01223_Pl_1")

#calculate bray-curtis distance 
funhel.dist=phyloseq::distance(sqrtfun, "bray") 

#remove this sample wich was low abundance in fungi
PS_normal=subset_samples(PS_normal, CollectionID !="WMEA_01223_Pl_1")

#sqrttransform
sqrt=transform_sample_counts(PS_normal, function(x) x^0.5)

#calculate bray curtis distance
hel.dist=phyloseq::distance(sqrt, "bray") 

#Calculate mantel
m=mantel(hel.dist,funhel.dist)
#make a dataframe for ggplot
mantel=as.data.frame(cbind(c(hel.dist), c(funhel.dist)))

# make plot
ggplot(mantel,aes(mantel$V1,mantel$V2))+
  geom_point()+
  geom_smooth()+
  theme_bw()+
  ggtitle("Fungi and Bacteria Community Similarity")+
  xlab("Fungal Dissimilarity")+
  ylab("Bacteria Dissimilarity")+
  annotate("text", x=.5, y=.3,label="r=0.434, P=0.001")+
  theme(plot.title = element_text(hjust = 1))

dev.print(png, "figures/community.mantel.png", units="in", width=4, height=4, res=300)
## quartz_off_screen 
##                 2
dev.off()
## null device 
##           1

Bacterial Abundance Occupancy

(problems here)

# Code by Anthony
# modified by Feresa Corazon, Chris Wall
# Purpose: Generate abundance occupancy graphs with bacteria and fungal data and site locations which are scaled the mean-annual precipitation 

# Load in packages
require('phyloseq')
require('bipartite')
require("raster")
library("ggplot2")
library("plotrix")
library("viridis")
library("lattice")


#Section 1A. Calling all data and values for bacteria
###################################
## Bacterial Abundance Occupancy ##
###################################

rarPhySeq=readRDS("data/rarPhySeq.gz")
#Calculate distance to shore
sample_data(rarPhySeq)$Shore_dist=pointDistance(cbind(sample_data(rarPhySeq)$Long,sample_data(rarPhySeq)$Lat), c(-158.062848, 21.640741),lonlat=TRUE)

# Merge by sample type
agg=merge_samples(rarPhySeq, "SampleType")

#Load in standarization code
phyloseq_standardize_otu_abundance <- function(physeq, method="total", ...){
  
  ## Check the orientation of the OTU table
  trows <- phyloseq::taxa_are_rows(physeq)
  if(trows == TRUE){ marg <- 2 } else { marg <- 1 }
  
  ## Extact OTU table
  comm <- as(object = phyloseq::otu_table(physeq), Class = "matrix")
  
  ## Standardize community table
  comm_std <- vegan::decostand(comm, method, MARGIN = marg, ...)
  
  ## Replace old otu_table with the new one
  phyloseq::otu_table(physeq) <- phyloseq::otu_table(comm_std, taxa_are_rows = trows)
  
  return(physeq)
}

#How many habitats is each ESV found? Standardize to convert to presence absence
habitats=colSums(otu_table(phyloseq_standardize_otu_abundance(agg, method = "pa")))
#habitats

#Merge by site location
sites=merge_samples(rarPhySeq, "FieldSite")

#Standardize to convert to presence absence and calculate site total
site=colSums(otu_table(phyloseq_standardize_otu_abundance(sites, method = "pa")))
#site

#Convert the otu table to presence absence in order to calculate range size
binarysite=phyloseq_standardize_otu_abundance(sites, method = "pa")

#Multiply each OTU by it's distance to shore
range=otu_table(binarysite)*sample_data(binarysite)$Shore_dist

#Convert 0 to NA
is.na(range) <- range==0

#Calculate min and max distance to shore
rangespan=apply(range,2,range, na.rm=TRUE)

#Subtract min distance from max distance
rangespan=rangespan[2,]-rangespan[1,]

#Do same thing but log transform range
nozerorangespan=rangespan
is.na(nozerorangespan) <- nozerorangespan==0

#Calculate total abundance
abund=colSums(otu_table(rarPhySeq))     

#Calculate how many samples an OTU is present
occupancy=colSums(otu_table(phyloseq_standardize_otu_abundance(rarPhySeq, method = "pa")))

#Calculate mean abundance per sample (where present)
meanabund=otu_table(rarPhySeq)
is.na(meanabund) <- meanabund==0
meanabund=colMeans(meanabund, na.rm=TRUE)

#Calculate per site habitat diversity (where present)
habpersite=cbind(colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s1"), method = "pa"))),
      colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s2"), method = "pa"))),
      colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s3"), method = "pa"))),
      colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s4"), method = "pa"))),
      colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s5"), method = "pa"))),
      colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s6"), method = "pa"))),
      colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s7"), method = "pa"))),
      colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s8"), method = "pa"))),
      colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s10"), method = "pa"))))
is.na(habpersite) <- habpersite==0
habpersite=rowMeans(habpersite, na.rm=TRUE)

Fungal Abundance Occupancy

# Code by Anthony
# modified by Feresa Corazon, Chris Wall
# Purpose: Generate abundance occupancy graphs with bacteria and fungal data and site locations which are scaled the mean-annual precipitation 


#Section 1B. Calling all data and values as above, but for fungi
################################
## Fungal Abundance Occupancy ##
################################
physeqF=readRDS("data/fungal_rarPhySeq.gz")
rarPhySeqF=rarefy_even_depth(physeqF, sample.size = 2000)

sample_data(rarPhySeqF)$Shore_dist=pointDistance(cbind(sample_data(rarPhySeqF)$Long,sample_data(rarPhySeqF)$Lat), c(-158.063625, 21.640741),lonlat=TRUE)
aggF=merge_samples(rarPhySeqF, "SampleType")
phyloseq_standardize_otu_abundance <- function(physeq, method="total", ...){
  
  ## Check the orientation of the OTU table
  trows <- phyloseq::taxa_are_rows(physeq)
  if(trows == TRUE){ marg <- 2 } else { marg <- 1 }
  
  ## Extact OTU table
  comm <- as(object = phyloseq::otu_table(physeq), Class = "matrix")
  
  ## Standardize community table
  comm_std <- vegan::decostand(comm, method, MARGIN = marg, ...)
  
  ## Replace old otu_table with the new one
  phyloseq::otu_table(physeq) <- phyloseq::otu_table(comm_std, taxa_are_rows = trows)
  
  return(physeq)
}

#How many habitats is each ESV found? Standardize to convert to presence absence
habitatsF=colSums(otu_table(phyloseq_standardize_otu_abundance(aggF, method = "pa")))
#habitatsF

#Merge by site location
sitesF=merge_samples(rarPhySeqF, "FieldSite")

#Standardize to convert to presence absence and calculate site total
siteF=colSums(otu_table(phyloseq_standardize_otu_abundance(sitesF, method = "pa")))
#siteF

#Convert the otu table to presence absence in order to calculate range size
binarysiteF=phyloseq_standardize_otu_abundance(sitesF, method = "pa")

#Multiply each OTU by it's distance to shore
rangeF=otu_table(binarysiteF)*sample_data(binarysiteF)$Shore_dist

#Convert 0 to NA
is.na(rangeF) <- rangeF==0

#Calculate min and max distance to shore
rangespanF=apply(rangeF,2,range.default, na.rm=TRUE)

#Subtract min distance from max distance
rangespanF=rangespanF[2,]-rangespanF[1,]

#Do same thing but log transform range
nozerorangespanF=rangespanF
is.na(nozerorangespanF) <- nozerorangespanF==0

#Calculate total abundance
abundF=colSums(otu_table(rarPhySeqF)) 

#Calculate how many samples an OTU is present
occupancyF=colSums(otu_table(phyloseq_standardize_otu_abundance(rarPhySeqF, method = "pa")))

#Calculate mean abundance per sample (where present)
meanabundF=otu_table(rarPhySeqF)
is.na(meanabundF) <- meanabundF==0
meanabundF=rowMeans(meanabundF, na.rm=TRUE)

#Calculate per site habitat diversity (where present)
habpersiteF=cbind(rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s1"), method = "pa"))),
                 rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s2"), method = "pa"))),
                 rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s3"), method = "pa"))),
                 rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s4"), method = "pa"))),
                 rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s5"), method = "pa"))),
                 rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s6"), method = "pa"))),
                 rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s7"), method = "pa"))),
                 rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s8"), method = "pa"))),
                 rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s10"), method = "pa"))))
is.na(habpersiteF) <- habpersiteF==0
habpersiteF=rowMeans(habpersiteF, na.rm=TRUE)
Bac-Fung: Habitat VS ESV range
#Section 2. Making all plots 
######################################################
########## ALL PLOTS AND STATISTICAL TESTS ###########
######################################################

# FIGURE 1. ESV RANGE RELATIONSHIPS
# Create a 3 x 3 plotting matrix
# The next  plots created will be plotted next to each other
par(mfrow = c(3, 2))


######  First row: plotting # OF HABITATS VS ESV RANGE
cor.test(rangespan, habitats) #Calculate correlation
## 
##  Pearson's product-moment correlation
## 
## data:  rangespan and habitats
## t = 74.414, df = 15724, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4986868 0.5218060
## sample estimates:
##       cor 
## 0.5103386
###### Plot Bacteria: # of Habitats vs. ESV Range (PLOT 1) ###### 
plot(jitter(rangespan,20),jitter(habitats,3),col="#481567FF", cex=.2, xlab="ESV range (m)", ylab="Number of Habitats ESV is Present", main="Bacteria: # of Habitats vs. ESV Range") + abline(lsfit(rangespan,habitats), col="black", lwd=3)
## integer(0)
ablineclip(lm(rangespan~habitats),lty = 2) #Add cut-off lines
text(4000,8, expression(paste(italic("cor = 0.512 \np <0.001"))), cex=0.8)


###### Plot Fungi: # of Habitats vs. ESV Range (PLOT 2) ###### 
cor.test(rangespanF, habitatsF) #Calculate correlation
## 
##  Pearson's product-moment correlation
## 
## data:  rangespanF and habitatsF
## t = 49.811, df = 3678, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6149939 0.6536001
## sample estimates:
##       cor 
## 0.6346928
plot(jitter(rangespanF,20),jitter(habitatsF,3),col="#95D840FF", cex=.2, xlab="ESV range (M)", ylab="Number of Habitats ESV is Present", main="Fungi: # of Habitats vs. ESV Range") + abline(lsfit(rangespanF,habitatsF), col="black", lwd=3)
## integer(0)
ablineclip(lm(rangespanF~habitatsF), lty = 2) #Add cut-off lines
text(4000,8, expression(paste(italic("cor = 0.625 \np <0.001"))), cex=0.8)


######  Second row: plotting MEAN ABUNDANCE VS ESV RANGE

######  Plot Bacteria: Mean Abundance vs. ESV Range (PLOT 3 ###### 
cor.test(rangespan, meanabund) #Calculate correlation
## 
##  Pearson's product-moment correlation
## 
## data:  rangespan and meanabund
## t = -0.46331, df = 15724, p-value = 0.6431
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.01932314  0.01193543
## sample estimates:
##          cor 
## -0.003694756
plot(jitter(rangespan,20),jitter(log(meanabund),3),col="#481567FF", cex=.2, xlab="ESV range (m)", ylab="Mean Abundance Per Sample (Where Present)", main="Bacteria: Mean Abundance vs. ESV Range")
abline(lsfit(rangespan,log(meanabund)), col="black", lwd=3)
ablineclip(lm(rangespan~meanabund), lty = 2) #Add cut-off lines
text(4000,8,expression(paste(italic("cor = -0.003 \np = 0.678"))), cex=0.8)

###### Plot Fungi: Mean abundance vs.  ESV Range (PLOT 4) ###### 
cor.test(rangespanF, meanabundF) #Calculate correlation
## 
##  Pearson's product-moment correlation
## 
## data:  rangespanF and meanabundF
## t = 5.3033, df = 3678, p-value = 1.204e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.05495683 0.11908890
## sample estimates:
##        cor 
## 0.08711312
plot(jitter(rangespanF,20),jitter(log(meanabundF),3),col="#95D840FF", cex=.2, xlab="ESV range (m)", ylab="Mean Abundance Per Sample (Where Present)", main=" Fungi: Mean Abundance vs. ESV Range")+ abline(lsfit(rangespanF,log(meanabundF)), col="black", lwd=3)
## integer(0)
ablineclip(lm(rangespanF~meanabundF), lty = 2) #Add cut-off lines
text(4000,8, expression(paste( italic("cor = 0.105 \np <0.001"))), cex=0.8)

###### Third row: plotting MEAN HABITAT OCCURENCE VS ESV RANGE ###### 

###### Plot Bacteria: Mean Habitat Occurence vs. ESV Range (PLOT 5) ###### 
cor.test(rangespan, habpersite) #Calculate correlation
## 
##  Pearson's product-moment correlation
## 
## data:  rangespan and habpersite
## t = 29.725, df = 15724, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2158068 0.2454031
## sample estimates:
##       cor 
## 0.2306583
plot(jitter(rangespan,20),jitter(habpersite,3),col="#481567FF", cex=.2, xlab="ESV range (m)", ylab="Mean Habitat Occurence (Per Site)", main="Bacteria: Mean Habitat Occurence vs. ESV Range") + abline(lsfit(rangespan,log(meanabund)), col="black", lwd=3) 
## integer(0)
ablineclip(lm(rangespan~habpersite), lty = 2) #Add cut-off lines
text(4000,8, expression(paste(italic("cor = 0.235 \np <0.001"))), cex=0.8)

###### Plot Fungi: Mean Habitat Occurence vs. ESV Range (PLOT 6) ###### 
cor.test(rangespanF, habpersiteF) #Calculate correlation
## 
##  Pearson's product-moment correlation
## 
## data:  rangespanF and habpersiteF
## t = 21.961, df = 3678, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3116030 0.3687402
## sample estimates:
##       cor 
## 0.3404859
plot(jitter(rangespanF,20),jitter(habpersiteF,3),col="#95D840FF", cex=.2, xlab="ESV range (m)", ylab="Mean Habitat Occurence (Per Site)", main="Fungi: Mean Habitat Occurence vs. ESV Range") + abline(lsfit(rangespanF,log(meanabundF)), col="black", lwd=3)
## integer(0)
ablineclip(lm(rangespanF~habpersiteF), lty = 2) #Add cut-off lines
text(4000,8, expression(paste( italic("cor = 0.333 \np <0.001"))))

dev.print(png, "figures/ESV.range.png", units="in", width=8, height=8, res=300)
## quartz_off_screen 
##                 2
dev.off()
## null device 
##           1

Rangespan of phytobiome communities

##FIGURE 2. RANGESPAN OF PHYTOBIOME COMMUNITIES 
#Set layout to combine boxplot and histogram
#Next plots will be created next to each other 
layout(mat = matrix(c(1,2,3,4),2,2),  height = c(1,8))

#First plot: Bacterial rangespan
#Set margin size for boxplot
par(mar=c(0, 4, 1.1, 2.1)) 
boxplot(rangespan, horizontal=TRUE, vertical=TRUE, xaxt="n" , col="#481567FF", main="Bacteria", frame=F) 

#Set margin size for histogram
par(mar=c(4, 4, 1.1, 2.1)) 
hist(rangespan, breaks=40 , col="#481567FF", main="", ylab = "Frequency", xlab="ESV Range (m)", ylim=c(0,5000), xlim=c(0,6000))

#Second Plot: Fungal rangespan
#Set margin size for boxplot
par(mar=c(0, 4, 1.1, 2.1))
boxplot(rangespanF, horizontal=TRUE , xaxt="n" , col="#95D840FF", main="Fungi", frame=F)

#Set margin size for histogram
par(mar=c(4, 4, 1.1, 2.1))
hist(rangespanF, breaks=40 , col="#95D840FF", main="", ylab = "Frequency", xlab="ESV Range (m)", ylim=c(0,3500), xlim=c(0,6000))

#Add title to figure!
title(main = "Rangespan of Phytobiome Communities", outer = TRUE, line = -0.5)

dev.print(png, "figures/community.range.png", units="in", width=6, height=6, res=300)
## quartz_off_screen 
##                 2
dev.off()
## null device 
##           1

Nestedness plots

#Code by Anthony
#To produce nestedness plots
require(phyloseq)
require(bipartite)
#read in rarefied Bacteria data
rarPhySeq=readRDS("data/rarPhySeq.gz")
#Merge samples by sample type
agg=merge_samples(rarPhySeq, "SampleType")
#convert to dataframe
aggtab=as.data.frame(otu_table(agg))
#calculate nested temperature
nested=nestedtemp(aggtab)

#read in Fungal data (not rarefied)
fungalphyseq1=readRDS("data/fungal_physeq1.gz")
#rarefy to 20000 sequences (drops 3 samples :())
funrarPhySeq=rarefy_even_depth(fungalphyseq1, sample.size=20000)
#merge by sample type
funagg=merge_samples(funrarPhySeq, "SampleType")
#convert to dataframe
funaggtab=as.data.frame(otu_table(funagg))
#calculate nested temp
funnested=nestedtemp(funaggtab)

#create EPS file
setEPS()
#using a postscript device
postscript("figures/nestedplots.eps")
#make a two row composite figure, reguce spacing between inner and outer margins
par(mfrow=c(2,1), mai=c(.1,1.,1,1), oma=c(.1,1,.1,1))
#plot the fungal nested figure, suppress taxon names
plot(funnested, kind="incid", names=c(TRUE,FALSE), col=c("white", "green"), lwd=3, main="Fungi Nested Plot")
#add stats as margin text
mtext("Nested Temp=42.5, P=0.001")
#plot the Bacteria nested figure, suppress taxon names
plot(nested, kind="incid", names=c(TRUE,FALSE), col=c("white", "purple"), lwd=3, main="Bacteria Nested Plot")
#add stats are margin text
mtext("Nested Temp=39.8, P=0.001")
dev.off()
## quartz_off_screen 
##                 2

Mantel Tests and Slope Plot

#######################################################################################
######################## Mantel Tests and Slope Plots By Plant Part ###################
################################ Chad Wilhite 04/25/19 ################################
#######################################################################################


# Load required packages
require(fossil)
require(vegan)
require(phyloseq)
require(ggplot2)


####################################
#### Mantel Tests by SampleType ####
####################################

#
##
###
#### Bacteria Analysis


# Load RDS file of normalized bacterial reads
PS_normal = readRDS("data/PS_normal.gz") #"PS_normal"


# Hellinger (square root) transform the data!
PS_normal.hell = transform_sample_counts(PS_normal, function(x) x^0.5)

# View SampleType names
unique(get_variable(PS_normal.hell, sample_variables(PS_normal.hell)[18]))
## [1] "Stem"    "Root"    "Air"     "Leaf"    "Soil"    "Litter"  "Axil"   
## [8] "Petiole"
#
##
### Subset out bacterial data by SampleType

# Subset bacterial data by plant part - the hard way... xD
stem_otu = subset_samples(PS_normal.hell, SampleType=="Stem")
root_otu = subset_samples(PS_normal.hell, SampleType=="Root")
air_otu  = subset_samples(PS_normal.hell, SampleType=="Air")
leaf_otu = subset_samples(PS_normal.hell, SampleType=="Leaf")
soil_otu = subset_samples(PS_normal.hell, SampleType=="Soil")
litt_otu = subset_samples(PS_normal.hell, SampleType=="Litter")
axil_otu = subset_samples(PS_normal.hell, SampleType=="Axil")
peti_otu = subset_samples(PS_normal.hell, SampleType=="Petiole")

#
##
### Generate Dissimilarity Matrices for bacterial OTU data

# Generate a dissimilarity matrix for grouped and individual plant part bacterial OTU data 
# adding "phyloseq::" to specify to run this through phyloseq and not other packages.

all_otu.dist  = phyloseq::distance(PS_normal.hell, "bray")
stem_otu.dist = phyloseq::distance(stem_otu, "bray") 
root_otu.dist = phyloseq::distance(root_otu, "bray") 
air_otu.dist  = phyloseq::distance(air_otu,  "bray") 
leaf_otu.dist = phyloseq::distance(leaf_otu, "bray") 
soil_otu.dist = phyloseq::distance(soil_otu, "bray")
litt_otu.dist = phyloseq::distance(litt_otu, "bray")
axil_otu.dist = phyloseq::distance(axil_otu, "bray") 
peti_otu.dist = phyloseq::distance(peti_otu, "bray")



#
##
### Create correctly ordered sized geographic distance dissimilarity matrix


#Gather bacterial spatial data from the phyloseq object
all.geo  = cbind(sample_data(PS_normal.hell)$Lon, sample_data(PS_normal.hell)$Lat)
stem.geo = cbind(sample_data(stem_otu)$Lon, sample_data(stem_otu)$Lat)
root.geo = cbind(sample_data(root_otu)$Lon, sample_data(root_otu)$Lat)
air.geo  = cbind(sample_data(air_otu)$Lon,  sample_data(air_otu)$Lat)
leaf.geo = cbind(sample_data(leaf_otu)$Lon, sample_data(leaf_otu)$Lat)
soil.geo = cbind(sample_data(soil_otu)$Lon, sample_data(soil_otu)$Lat)
litt.geo = cbind(sample_data(litt_otu)$Lon, sample_data(litt_otu)$Lat)
axil.geo = cbind(sample_data(axil_otu)$Lon, sample_data(axil_otu)$Lat) 
peti.geo = cbind(sample_data(peti_otu)$Lon, sample_data(peti_otu)$Lat)



# Generate a dissimilarity matrix for grouped (all plant parts) and by plant part 
#   bacterial geographic distance data 
all.geodist  = earth.dist(all.geo)
stem.geodist = earth.dist(stem.geo)
root.geodist = earth.dist(root.geo)
air.geodist  = earth.dist(air.geo)
leaf.geodist = earth.dist(leaf.geo)
soil.geodist = earth.dist(soil.geo)
litt.geodist = earth.dist(litt.geo)
axil.geodist = earth.dist(axil.geo)
peti.geodist = earth.dist(peti.geo)



#Mantel test by SampleType
all_bacteria_mantel  = mantel(log( all.geodist+1),log( all_otu.dist), permutations= 999)
stem_bacteria_mantel = mantel(log(stem.geodist+1),log(stem_otu.dist), permutations= 999)
root_bacteria_mantel = mantel(log(root.geodist+1),log(root_otu.dist), permutations= 999)
air_bacteria_mantel  = mantel(log( air.geodist+1),log( air_otu.dist), permutations= 999)
leaf_bacteria_mantel = mantel(log(leaf.geodist+1),log(leaf_otu.dist), permutations= 999)
soil_bacteria_mantel = mantel(log(soil.geodist+1),log(soil_otu.dist), permutations= 999)
litt_bacteria_mantel = mantel(log(litt.geodist+1),log(litt_otu.dist), permutations= 999)
axil_bacteria_mantel = mantel(log(axil.geodist+1),log(axil_otu.dist), permutations= 999)
peti_bacteria_mantel = mantel(log(peti.geodist+1),log(peti_otu.dist), permutations= 999)




#Call the bacterial results
all_bacteria_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(all.geodist + 1), ydis = log(all_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.1211 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0427 0.0560 0.0678 0.0878 
## Permutation: free
## Number of permutations: 999
stem_bacteria_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(stem.geodist + 1), ydis = log(stem_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.2521 
##       Significance: 0.184 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.322 0.371 0.400 0.443 
## Permutation: free
## Number of permutations: 999
root_bacteria_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(root.geodist + 1), ydis = log(root_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.445 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.150 0.209 0.276 0.329 
## Permutation: free
## Number of permutations: 999
air_bacteria_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(air.geodist + 1), ydis = log(air_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.2071 
##       Significance: 0.256 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.326 0.371 0.421 0.463 
## Permutation: free
## Number of permutations: 999
leaf_bacteria_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(leaf.geodist + 1), ydis = log(leaf_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.2337 
##       Significance: 0.151 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.409 0.505 0.567 0.626 
## Permutation: free
## Number of permutations: 999
soil_bacteria_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(soil.geodist + 1), ydis = log(soil_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.04483 
##       Significance: 0.382 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.397 0.496 0.574 0.630 
## Permutation: free
## Number of permutations: 999
litt_bacteria_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(litt.geodist + 1), ydis = log(litt_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.2802 
##       Significance: 0.135 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.317 0.386 0.441 0.503 
## Permutation: free
## Number of permutations: 999
axil_bacteria_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(axil.geodist + 1), ydis = log(axil_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.07522 
##       Significance: 0.375 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.367 0.454 0.516 0.573 
## Permutation: free
## Number of permutations: 999
peti_bacteria_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(peti.geodist + 1), ydis = log(peti_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.2122 
##       Significance: 0.165 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.381 0.511 0.567 0.619 
## Permutation: free
## Number of permutations: 999
#############################################
################ Fungal #####################
#############################################

#Load RDS file of fungal reads
FPS_normal = readRDS("data/fungal_PS_normal.gz") #"fungal_PS_normal"

#square root transform (Hellinger)
FPS_normal.hell = transform_sample_counts(FPS_normal, function(x) x^0.5)



#Subset Fungal Data by SampleType
fstem_otu = subset_samples(FPS_normal.hell, SampleType=="Stem")
froot_otu = subset_samples(FPS_normal.hell, SampleType=="Root")
fair_otu  = subset_samples(FPS_normal.hell, SampleType=="Air")
fleaf_otu = subset_samples(FPS_normal.hell, SampleType=="Leaf")
fsoil_otu = subset_samples(FPS_normal.hell, SampleType=="Soil")
flitt_otu = subset_samples(FPS_normal.hell, SampleType=="Litter")
faxil_otu = subset_samples(FPS_normal.hell, SampleType=="Axil")
fpeti_otu = subset_samples(FPS_normal.hell, SampleType=="Petiole")


#Generate Distance Matrix
fall_otu.dist  = phyloseq::distance(FPS_normal.hell, "bray")
fstem_otu.dist = phyloseq::distance(fstem_otu, "bray")
froot_otu.dist = phyloseq::distance(froot_otu, "bray")
fair_otu.dist  = phyloseq::distance( fair_otu, "bray")
fleaf_otu.dist = phyloseq::distance(fleaf_otu, "bray")
fsoil_otu.dist = phyloseq::distance(fsoil_otu, "bray")
flitt_otu.dist = phyloseq::distance(flitt_otu, "bray")
faxil_otu.dist = phyloseq::distance(faxil_otu, "bray")
fpeti_otu.dist = phyloseq::distance(fpeti_otu, "bray")




#Gather fungal spatial data from the phyloseq object
fall.geo  = cbind(sample_data(FPS_normal.hell)$Lon, sample_data(FPS_normal.hell)$Lat)
fstem.geo = cbind(sample_data(fstem_otu)$Lon, sample_data(fstem_otu)$Lat)
froot.geo = cbind(sample_data(froot_otu)$Lon, sample_data(froot_otu)$Lat)
fair.geo  = cbind(sample_data( fair_otu)$Lon, sample_data( fair_otu)$Lat)
fleaf.geo = cbind(sample_data(fleaf_otu)$Lon, sample_data(fleaf_otu)$Lat)
fsoil.geo = cbind(sample_data(fsoil_otu)$Lon, sample_data(fsoil_otu)$Lat)
flitt.geo = cbind(sample_data(flitt_otu)$Lon, sample_data(flitt_otu)$Lat)
faxil.geo = cbind(sample_data(faxil_otu)$Lon, sample_data(faxil_otu)$Lat)
fpeti.geo = cbind(sample_data(fpeti_otu)$Lon, sample_data(fpeti_otu)$Lat)


#Create correct sized fungal geographic distance dissimilarity matrix
fall.geodist  = earth.dist( fall.geo)
fstem.geodist = earth.dist(fstem.geo)
froot.geodist = earth.dist(froot.geo)
fair.geodist  = earth.dist( fair.geo)
fleaf.geodist = earth.dist(fleaf.geo)
fsoil.geodist = earth.dist(fsoil.geo)
flitt.geodist = earth.dist(flitt.geo)
faxil.geodist = earth.dist(faxil.geo)
fpeti.geodist = earth.dist(fpeti.geo)





#Mantel test by SampleType
fall_mantel  = mantel(log( fall.geodist+1),log( fall_otu.dist), permutations= 999)
fstem_mantel = mantel(log(fstem.geodist+1),log(fstem_otu.dist), permutations= 999)
froot_mantel = mantel(log(froot.geodist+1),log(froot_otu.dist), permutations= 999)
fair_mantel  = mantel(log( fair.geodist+1),log( fair_otu.dist), permutations= 999)
fleaf_mantel = mantel(log(fleaf.geodist+1),log(fleaf_otu.dist), permutations= 999)
fsoil_mantel = mantel(log(fsoil.geodist+1),log(fsoil_otu.dist), permutations= 999)
flitt_mantel = mantel(log(flitt.geodist+1),log(flitt_otu.dist), permutations= 999)
faxil_mantel = mantel(log(faxil.geodist+1),log(faxil_otu.dist), permutations= 999)
fpeti_mantel = mantel(log(fpeti.geodist+1),log(fpeti_otu.dist), permutations= 999)




#call the fungal results
fall_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(fall.geodist + 1), ydis = log(fall_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.3641 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0691 0.0897 0.1127 0.1427 
## Permutation: free
## Number of permutations: 999
fstem_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(fstem.geodist + 1), ydis = log(fstem_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.5552 
##       Significance: 0.002 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.278 0.339 0.382 0.426 
## Permutation: free
## Number of permutations: 999
froot_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(froot.geodist + 1), ydis = log(froot_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.6999 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.274 0.336 0.371 0.438 
## Permutation: free
## Number of permutations: 999
fair_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(fair.geodist + 1), ydis = log(fair_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.6131 
##       Significance: 0.003 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.339 0.421 0.484 0.540 
## Permutation: free
## Number of permutations: 999
fleaf_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(fleaf.geodist + 1), ydis = log(fleaf_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.1138 
##       Significance: 0.235 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.239 0.351 0.397 0.476 
## Permutation: free
## Number of permutations: 999
fsoil_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(fsoil.geodist + 1), ydis = log(fsoil_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.6737 
##       Significance: 0.006 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.418 0.509 0.583 0.645 
## Permutation: free
## Number of permutations: 999
flitt_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(flitt.geodist + 1), ydis = log(flitt_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.5981 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.336 0.419 0.457 0.518 
## Permutation: free
## Number of permutations: 999
faxil_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(faxil.geodist + 1), ydis = log(faxil_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.2254 
##       Significance: 0.165 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.284 0.353 0.398 0.436 
## Permutation: free
## Number of permutations: 999
fpeti_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = log(fpeti.geodist + 1), ydis = log(fpeti_otu.dist),      permutations = 999) 
## 
## Mantel statistic r: 0.6702 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.294 0.376 0.436 0.514 
## Permutation: free
## Number of permutations: 999
###################################################
# Save Mantel Correlation Statistics and P Values #
###################################################

#Call plant part results
mant_cor = c(stem_bacteria_mantel$statistic, root_bacteria_mantel$statistic, 
    air_bacteria_mantel$statistic, leaf_bacteria_mantel$statistic, 
    soil_bacteria_mantel$statistic, litt_bacteria_mantel$statistic, 
    axil_bacteria_mantel$statistic, peti_bacteria_mantel$statistic, 
    fstem_mantel$statistic, froot_mantel$statistic, fair_mantel$statistic, 
    fleaf_mantel$statistic, fsoil_mantel$statistic, flitt_mantel$statistic, 
    faxil_mantel$statistic, fpeti_mantel$statistic)

#Call plant part p-values
p_val = c(stem_bacteria_mantel$signif, root_bacteria_mantel$signif, 
    air_bacteria_mantel$signif, leaf_bacteria_mantel$signif, 
    soil_bacteria_mantel$signif, litt_bacteria_mantel$signif, 
    axil_bacteria_mantel$signif, peti_bacteria_mantel$signif, 
    fstem_mantel$signif, froot_mantel$signif, fair_mantel$signif, 
    fleaf_mantel$signif, fsoil_mantel$signif, flitt_mantel$signif, 
    faxil_mantel$signif, fpeti_mantel$signif)

#Create results data frame with mantel R and p-values for each plant part
results = data.frame( c( rep('Bacterial',8), rep('Fungal',8) ), 
    c( rep( c('Stem', 'Root', 'Air', 'Leaf', 'Soil', 'Litter', 'Axil', 'Petiole'),2 ) ),
    mant_cor, p_val
    )

#Name the columns intelligible names
colnames(results) = c('type', 'Part', 'mant_cor', 'p_val')

#Adjust p-values for repeated tests on the within plant part groups for bacteria
bact.cor.p = p.adjust(results[results$type == 'Bacterial',4], method = 'bonferroni')
bact.cor.p = data.frame(cor.p = bact.cor.p, type = rep('Bacterial',8), 
    Part = as.character(unique(results$Part)) )

#Adjust p-values for repeated tests on the within plant part groups for fungi
fung.cor.p = p.adjust(results[results$type == 'Fungal',4], method = 'bonferroni')
fung.cor.p = data.frame(cor.p = fung.cor.p, type = rep('Fungal',8), 
    Part = as.character(unique(results$Part)) )

#Join corrected p-values to for fungus and bacteria together
cor.p = rbind(bact.cor.p, fung.cor.p)


#Join corrected p-values to results data frame
results = merge(results, cor.p)



#Add overall test (all parts together) to data frame
results2 = rbind( results, data.frame( type = c('Bacterial', 'Fungal'), Part = rep('All',2),
    mant_cor = c(all_bacteria_mantel$statistic, fall_mantel$statistic), 
    p_val = c(all_bacteria_mantel$signif, fall_mantel$signif),
    cor.p = rep('NA', 2))
    )

#Save all result data as a csv
write.csv(results2, file = 'output/Mantel_Plant_Part_Results_Hell_Bray.csv')




####################################################
################### Slopes plot ####################
####################################################


# The colorblind friendly palette with grey:
cbPalette = c('#999999', '#E69F00', '#56B4E9', 
    '#009E73', '#F0E442', '#0072B2', '#D55E00', '#CC79A7')



###################
#### Bacterial ####
###################

#Combine all plant part bacterial OTU dissimilarity matrices
ball.otudist = c(stem_otu.dist, root_otu.dist, air_otu.dist,
    leaf_otu.dist, soil_otu.dist, litt_otu.dist, axil_otu.dist, peti_otu.dist)


#Combine all plant part bacterial geographic dissimilarity matrices
ball.geodist = c(stem.geodist, root.geodist, air.geodist,
    leaf.geodist, soil.geodist, litt.geodist, axil.geodist, peti.geodist)

#Set up properly ordered names by plant part
bplpt = c( rep('Stem', 36), rep('Root', 36), rep('Air', 36), rep('Leaf', 36),
     rep('Soil', 36), rep('Litter', 36), rep('Axil', 36), rep('Petiole', 36) ) 


#Combine the OTU and geographic dissimilarity vectors with their plant part names
bac.all.dist = data.frame(b.geo.dist = ball.geodist, b.otu.dist = ball.otudist, b.pl.part = bplpt)


#Give each plant part our particular colorblind friendly color
b.part.col = data.frame(cbPalette, levels(bac.all.dist$b.pl.part))
colnames(b.part.col) = c('b.col', 'b.pl.part')
b.part.col$b.col = as.character(b.part.col$b.col)

#Merge assigned plant part colors in the main data frame
bac.all.dist = merge(bac.all.dist, b.part.col)


#Create and save a plot of slopes for each plant part!
bac.slope.plot = ggplot(bac.all.dist, aes(x=b.geo.dist, y=b.otu.dist, color = b.pl.part)) +
    scale_color_manual( values = unique(bac.all.dist$b.col), name = "Sample Type") +
    geom_smooth(method = lm, se = FALSE) + geom_point() +
    theme(axis.line = element_line(color = "black"), legend.background = element_rect(),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), legend.key=element_blank(),
        legend.justification=c(1,0), legend.position=c(1,0)) +
    scale_x_continuous(name = "Pairwise Geographic Dissimilarity", breaks = 0:5) + 
    ylab("Pairwise Community Dissimilarity")

#Show plot
bac.slope.plot

#Save plot
ggsave(filename = "figures/Bacterial_Slopes_by_Plant_Part.eps", width = 10, height = 10)

dev.off()
## null device 
##           1
################
#### Fungal ####
################

#Combine all fungal plant part OTU dissimilarity matrices
fall.otudist = c(fstem_otu.dist, froot_otu.dist, fair_otu.dist,
    fleaf_otu.dist, fsoil_otu.dist, flitt_otu.dist, faxil_otu.dist, fpeti_otu.dist)


#Combine all fungal plant part geographic dissimilarity matrices
fall.geodist = c(fstem.geodist, froot.geodist, fair.geodist,
    fleaf.geodist, fsoil.geodist, flitt.geodist, faxil.geodist, fpeti.geodist)

#Set up properly ordered names by plant part
plpt = c( rep('Stem', 36), rep('Root', 36), rep('Air', 36), rep('Leaf', 36),
     rep('Soil', 36), rep('Litter', 36), rep('Axil', 36), rep('Petiole', 36) ) 


#Combine the fungal OTU and geographic dissimilarity vectors with their plant part names
fun.all.dist = data.frame(geo.dist = fall.geodist, otu.dist = fall.otudist, pl.part = plpt)


#Give each plant part our particular colorblind friendly color
part.col = data.frame(cbPalette, levels(fun.all.dist$pl.part))
colnames(part.col) = c('col', 'pl.part')
part.col$col = as.character(part.col$col)

fun.all.dist = merge(fun.all.dist, part.col)


#Make and save a plot of slopes for each fungal plant part!
fun.slope.plot = ggplot(fun.all.dist, aes(x=geo.dist, y=otu.dist, color = pl.part)) +
    scale_color_manual( values = unique(fun.all.dist$col), name = "Sample Type") +
    geom_smooth(method = lm, se = FALSE) + geom_point() +
    theme(axis.line = element_line(color = "black"), legend.background = element_rect(),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), legend.key=element_blank(),
        legend.justification=c(1,0), legend.position=c(1,0)) +
     scale_x_continuous(name = "Pairwise Geographic Dissimilarity", breaks = 0:5) + 
    ylab("Pairwise Community Dissimilarity")

#Show plot
fun.slope.plot

#Save plot
ggsave(filename = "figures/Fungal_Slopes_by_Plant_Part.eps", width = 10, height = 10)

dev.off()
## null device 
##           1